#Jaime E. Settle
#1.27.16
#Political Chameleons Study 1 Replication Code

#This file includes the replication code for analyses in Study 1 of "Political Chameleons: An Exploration of Conformity in Political Discussions"
#The file first shows the analyses for summary statistics, manipulation checks, and treatment main effects. Second, it includes the code that generated the figures in the paper. Third, it includes the code for the regression models

rm(list=ls())

#Load packages 
library(stargazer)
library(Zelig)
library(ZeligChoice)
library(xtable)

#Read in Data
study1data <-read.csv("study1data.csv")
dim(study1data)

#########################################################################
#########################################################################
###------------------Results Presented in Main Text-------------------###
#########################################################################
#########################################################################


#--------Summary Statistics-------------#
#See the code for the Appendix for the figures and balance tables referenced in the text

#What is the likelihood that Sally expresses her true opinion to the group? (value of 1 or 2)
table(study1data$Truefeelings)

P1 <- length(which(study1data$Truefeelings<3))
N1 <- length(which(!is.na(study1data$Truefeelings)))
SAMPLEPROP1 <- P1/N1
SQRT1 <- sqrt((SAMPLEPROP1*(1-SAMPLEPROP1))/N1)
SEU1 <- SAMPLEPROP1 + 1.96*SQRT1
SEL1 <- SAMPLEPROP1 - 1.96*SQRT1

#Social Consequence Variables Distribution (value of 3 or lower)
P1 <- length(which(study1data$NewJob<4))
N1 <- length(which(!is.na(study1data$NewJob)))
SAMPLEPROP1 <- P1/N1

P1 <- length(which(study1data$Dinner<4))
N1 <- length(which(!is.na(study1data$Dinner)))
SAMPLEPROP1 <- P1/N1


#--------T.Test Results-------------#

#Comparing expression of true opinion in Partisan Minority vs. Balanced Party
t.test(study1data$Truefeelings ~ study1data$PartisanMinority)

#Comparing true opinion proportions between conditions
prop.test(c(length(which(study1data$Truefeelings==1 & study1data$PartisanMinority==1)), length(which(study1data$Truefeelings==1 & study1data$PartisanMinority==0))), c(length(which(study1data$PartisanMinority==1)), length(which(study1data$PartisanMinority==0))))

#Comparing New Job outcomes
t.test(study1data$NewJob ~ study1data$PartisanMinority)

#Comparing New Job proportions between conditions
prop.test(c(length(which(study1data$NewJob>2 & study1data$PartisanMinority==1)), length(which(study1data$NewJob>2 & study1data$PartisanMinority==0))), c(length(which(study1data$PartisanMinority==1)), length(which(study1data$PartisanMinority==0))))

#-------- Figure 1-------------#

# Uncomfortable
c <- t.test(study1data$Uncomfortable ~ study1data$PartisanMinority)
d <- t.test(study1data$Uncomfortable ~ study1data$OpposeSally)

# True Feelings
a <- t.test(study1data$Truefeelings ~ study1data$PartisanMinority)
b <- t.test(study1data$Truefeelings ~ study1data$OpposeSally)

#Dinner
e <- t.test(study1data$Dinner ~ study1data$PartisanMinority)
f <- t.test(study1data$Dinner ~ study1data$OpposeSally)

# NewJob
g <- t.test(study1data$NewJob ~ study1data$PartisanMinority)
h <- t.test(study1data$NewJob ~ study1data$OpposeSally)

m <- matrix(NA,8,3)
m[1,1] <- (a$estimate[2] - a$estimate[1])
m[2,1] <- (b$estimate[2] - b$estimate[1])
m[3,1] <- (c$estimate[2] - c$estimate[1])
m[4,1] <- (d$estimate[2] - d$estimate[1])
m[5,1] <- (e$estimate[2] - e$estimate[1])
m[6,1] <- (f$estimate[2] - f$estimate[1])
m[7,1] <- (g$estimate[2] - g$estimate[1])
m[8,1] <- (h$estimate[2] - h$estimate[1])

m[1,2] <- a$conf.int[1]*-1
m[2,2] <- b$conf.int[1]*-1
m[3,2] <- c$conf.int[1]*-1
m[4,2] <- d$conf.int[1]*-1
m[5,2] <- e$conf.int[1]*-1
m[6,2] <- f$conf.int[1]*-1
m[7,2] <- g$conf.int[1]*-1
m[8,2] <- h$conf.int[1]*-1

m[1,3] <- a$conf.int[2]*-1
m[2,3] <- b$conf.int[2]*-1
m[3,3] <- c$conf.int[2]*-1
m[4,3] <- d$conf.int[2]*-1
m[5,3] <- e$conf.int[2]*-1
m[6,3] <- f$conf.int[2]*-1
m[7,3] <- g$conf.int[2]*-1
m[8,3] <- h$conf.int[2]*-1

m

seq <- 8:1
seq2 <- c(1.5, 3.5, 5.5, 7.5)
labeltext <- c( "Likelihood that Sally\nwants to look for\na new job within\nthe next year?", "Likelihood that Sally\ninvites her coworkers over\nfor dinner in the\nnext six months?", "Likelihood that Sally\nfeels uncomfortable\nanswering this question?", "Likelihood that Sally\nexpresses her true\nopinion to the group?")
labelcolor <- c("gray55", "gray55", "gray55", "gray 20")
color <- rep(c("chocolate2", "orchid3"), 4)
pchcode <- rep(c(15, 17), 4)

par(mar=c(3,15,3,1), col=c("gray40"), col.axis=c("gray40"), col.lab=c("gray40"), col.main=c("black"))
plot(m[,1], seq, type="n", col=color, pch=pchcode, xlim=c(-.5, .5), cex=1.5, yaxt="n", xaxt="n", ylab="", xlab="", font=2, font.lab=2, cex.lab=1.5, main="Mean Differences\nbetween Treatment Conditions")

for(i in 1:8){
lines(c(m[i,2], m[i,3]), c(seq[i], seq[i]), col=c("grey60"))
points(m[i,1], seq[i], col=color[i], pch=pchcode[i], cex=1.5, )
}
mtext(side=2, at=c(seq2), text=labeltext, font=2, cex=1.25, las=1, col=labelcolor, line = 1)
axis(side=1, at=c(-0.5, -.25, 0, .25, .5), labels=c("-0.50", "-0.25", "0", "0.25", "0.50"), tick=T, col.ticks=c("gray40"), col.lab=c("black"), col.axis=c("black"))

legend(-.5,2, legend=c("Contextual Manipulation", "Provocation Manipulation"), cex=.7, col=color, pch=pchcode)
abline(v=0, lty=2, col="grey60")
abline(h=6.5, lty=1, col="grey20")


####------------------Regression Results and Table 1---------------####

model1 <- zelig(as.factor(Truefeelings) ~ Uncomfortable, model="ologit", data=study1data)
summary(model1)
ctable <- coef(summary(model1))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable1 <- cbind(ctable, `p value` = p)
ctable1

model2 <- zelig(as.factor(Truefeelings) ~ Uncomfortable + male + white, model="ologit", data=study1data)
summary(model2)
ctable <- coef(summary(model2))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable2 <- cbind(ctable, `p value` = p)
ctable2

model3 <- zelig(as.factor(Truefeelings) ~ Uncomfortable + male + white + engage2012 +polknow + polinterest_R + pidstrength + ideo7, model="ologit", data=study1data)
summary(model3)
ctable <- coef(summary(model3))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable3 <- cbind(ctable, `p value` = p)
ctable3


stargazer(ctable1, ctable2, ctable3)
stargazer(model1, model2, model3)



#########################################################################
#########################################################################
###------------------------Appendix Figures----------------------------###
#########################################################################
#########################################################################

###---------------Table 5 in the Appendix--------------------###

#Balance Table for Fall vs. Spring Data Collection
balvar <- c(10,9,20,14,18,15,16,19,12)
baltable <- matrix(NA, length(balvar), 4)

for (i in 1:length(balvar)){

baltable[i,1] <- names(study1data)[balvar[i]]
baltable[i,2] <- round(mean(study1data[which(study1data$semester=="f13"), balvar[i]], na.rm=T), 3)
baltable[i,3] <- round(mean(study1data[which(study1data$semester=="S14"), balvar[i]], na.rm=T), 3)
baltest <- t.test(study1data[,balvar[i]] ~ study1data$semester)
baltable[i,4] <-round(baltest$p.value, 3)

}

table(study1data$semester)

xtable(baltable, digits=3)

#Balance Table for Treatment Groups
balvar <- c(10,9,20,14,18,15,16,19,12)
baltable <- matrix(NA, length(balvar), 4)

for (i in 1:length(balvar)){

baltable[i,1] <- names(study1data)[balvar[i]]
baltable[i,2] <- round(mean(study1data[which(study1data$PartisanMinority==1), balvar[i]], na.rm=T), 3)
baltable[i,3] <- round(mean(study1data[which(study1data$PartisanMinority==0), balvar[i]], na.rm=T), 3)
baltest <- t.test(study1data[,balvar[i]] ~ study1data$PartisanMinority)
baltable[i,4] <-round(baltest$p.value, 3)

}

table(study1data$PartisanMinority)

xtable(baltable, digits=3)

###--------------Figure 6 in the Appendix--------------------###

###Histogram for True Feelings###

barplot(table(study1data$Truefeelings), space=0,
     beside=T,col=c("blue"),font=2,font.lab=2,ylim=c(0,200), ylab="", yaxt="n",
    
)

axis(side=2, at=seq(from=0, to=200, by=50), label=c("0", "50", "100", "150", "200"), las=2, cex=.5)
axis(side=1, at=c(2.5), label=c("Likelihood of Reporting True Feelings"), line=1.5, lwd=0, font=2)
axis(side=1, at=c(2.5), label=c("5='Very Likely'"), line=3, lwd=0)

axis(side=2, at=100, label=c("Count"), las=0, cex=.7, line=2, font=2, tick=F)
mtext(c("Likelihood of Expressing True Feelings"), las=0, cex=1.5, line=1, font=2, side=3, at=2.5)

location <- c(.5, 1.5, 2.5, 3.5, 4.5)
cvalue <- table(study1data$Truefeelings)/length(!is.na(study1data$Truefeelings))
yloc <- c(50, 125, 165, 125, 50)

for (i in 1:length(cvalue)){
text(x=(location[i]), y=yloc[i], labels=paste(round(cvalue[i]*100, 2), "%", sep=""), font=2, cex=1.1, col=c("gray30"))
}

abline(h=0, lty=3, col=c("gray30"))

###Histogram for Uncomfortabe###
barplot(table(study1data$Uncomfortable), space=0,
     beside=T,col=c("seagreen3"),font=2,font.lab=2,ylim=c(0,200), ylab="", yaxt="n",
    
)

axis(side=2, at=seq(from=0, to=200, by=50), label=c("0", "50", "100", "150", "200"), las=2, cex=.5)
axis(side=1, at=c(2.5), label=c("Likelihood of Feeling Uncomfortable"), line=1.5, lwd=0, font=2)
axis(side=1, at=c(2.5), label=c("5='Very Likely'"), line=3, lwd=0)

axis(side=2, at=100, label=c("Count"), las=0, cex=.7, line=2, font=2, tick=F)
mtext(c("Likelihood of Feeling Uncomfortable"), las=0, cex=1.5, line=1, font=2, side=3, at=2.5)

location <- c(.5, 1.5, 2.5, 3.5, 4.5)
cvalue <- table(study1data$Uncomfortable)/length(!is.na(study1data$Uncomfortable))
yloc <- c(25, 40, 85, 165, 175)

for (i in 1:length(cvalue)){
text(x=(location[i]), y=yloc[i], labels=paste(round(cvalue[i]*100, 2), "%", sep=""), font=2, cex=1.1, col=c("gray30"))
}

abline(h=0, lty=3, col=c("gray30"))


####-----------------Table 6 in the Appendix---------------####

testA <- t.test(study1data$Uncomfortable ~ study1data$Truefeelings < 3)
testB <- t.test(study1data$NewJob ~ study1data$Truefeelings < 3)
testC <- t.test(study1data$Dinner ~ study1data$Truefeelings < 3)

sumstat <- matrix(NA, 3,3)

sumstat[1,1] <- testA$estimate[1]
sumstat[1,2] <- testA$estimate[2]
sumstat[1,3] <- testA$p.value
sumstat[2,1] <- testB$estimate[1]
sumstat[2,2] <- testB$estimate[2]
sumstat[2,3] <- testB$p.value
sumstat[3,1] <- testC$estimate[1]
sumstat[3,2] <- testC$estimate[2]
sumstat[3,3] <- testC$p.value

xtable(sumstat, digits=3)
